tic
% Построение реализаций при гармоническом воздействии
% 	Фролов Д.А. ЭР-23-09
clc       % очистка экрана
close all % закрытие всех лишних окон
L=100e-6; % Катушка индуктивности, Гн
global L  % Глобальная переменная, чтобы использовать в функции f
R=47;	  % Сопротивление резистора в цепи, Ом
global R  % Глобальная переменная, чтобы использовать в функции f
A=1;      % амплитуда входного гармонического воздействия, В
global A  % Глобальная переменная, чтобы использовать в функции f
w0=0;     % Частота анализа
global w0 % Глобальная переменная, чтобы использовать в функции f

% Для решения ДУ в мат. пакете Octave есть функция lsode(f,x0,t)
% f  - имя функции, которая описывает ДУ, имеет два входных параметра
% х0 - начальное условие в первый момент времени
% t  - вектор, который содержит отсчеты времени. Первый элемент - для начального условия

function xdot=f(x,t)
	global w0 % используем заранее объявленную частоту
	global L  % используем заранее объявленную индуктивность
	global R  % используем заранее объявленное сопротивление
	global A  % используем заранее объявленную амплитуду
	xdot=-2*R/L*x+A*cos(w0*t)/L; % pi=-2R/L * i + uвх(t)
endfunction

Nt=10000;          % количество отсчетов по времени
tmin=0;            % начало временного интервала анализа
tmax=75e-6;        % конец временного анализа, 75мкс, если взять 50, то на малых частотах получим скачок
dt=(tmax-tmin)/Nt; % шаг по времени, 7.5 нс
t=(0:Nt-1)*dt;     % задаю вектор времени

Nw=5;              % количество отсчетов по частоте
wmin=0;            % начало частотного интервала анализа
wmax=10e6;         % конец частотного анализа
dw=(wmax-wmin)/Nw; % шаг по частоте, 2 МГц
w=(0:Nw-1)*dw;     % задаю вектор частоты

x0=0;              % задаю начальное условие в нулевой момент времени

for i=1:Nw                  % цикл от 5 раз
	w0=w(i);                % перебираем все значения частоты
	cur=lsode('f',x0,t);    % рещаем ДУ с новым значение w0
	figure(i)               % для i-ой плоскости построения графиков
							% построить реализацию входного воздействия и выходного
	plot(t/1e-6,A*cos(w0*t),t/1e-6,cur*R)
	xlabel('t, us')         % подпись оси ОХ
	ylabel('Uin, Uout, V')  % подпись оси OY
	title('OSC')     		% Название графика
	grid on                 % вывести сетку
	box off                 % оставить только оси
	legend('Uin','Uout')    % вывести легенду
	ylim([-1,1])            % границы отображения по оси OY
	if (w0==0)              % если частота нулевая 
		xlim([0,10])        % то по времени границы ставим от 0 до 10 мкс
		else                % иначе по времени границы ставим от 0 до 2ух периодов входного воздействия
		xlim([0,2*pi/abs(w0)*2/1e-6]) % гпаницы отображения по оси ОХ
		endif
	end
toc